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Abstract 

Stable  and  spectrally  accurate  numerical  methods  are  constructed  on  arbitrary  grids 
for  partial  differential  equations.  These  new  methods  are  equivalent  to  conventional  spectral 
methods  but  do  not  rely  on  specific  grid  distributions.  Specifically,  we  show  how  to  implement 
Legendre  Galerkin,  Legendre  collocation,  and  Laguerre  Galerkin  methodology  on  arbitrary 
grids. 
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1  Introduction 


Conventional  spectral  methods  impose  rigid  requirements  on  the  computational  grids  used. 
The  grid  points  are  the  nodes  of  Gauss-like  quadrature  formulas  (Gauss,  Gauss  Radau,  or 
Gauss  Lobatto  (GL)  formulas).  These  nodes  are  denser  at  the  boundaries  than  in  the  middle 
of  the  domain.  Although  this  property  is  suitable  for  boundary-layer  problems,  it  may  create 
difficulties  for  other  types  of  problems,  particularly  those  with  disparate  length  scales  that 
occur  in  multiple  regions  of  the  domain  (e.g.,  diffusive  burning  or  detonation  and  reacting 
mixing  layers).  The  principle  reason  for  the  degradation  in  performance  on  these  disparate 
problems  is  that  the  predetermined  node  points  do  not,  in  general,  coincide  with  the  features 
that  are  being  resolved.  Extensive  mappings  can  concentrate  the  node  points  into  regions 
more  ideally  suited  for  accurate  resolution  but  present  a  serious  limitation  for  complicated 
problems.  For  this  reason,  spectral  multidomain  techniques  have  an  obvious  advantage  for 
complicated  problems  [l]-[3]. 

Another  complication  that  conventional  spectral  methods  have,  is  their  implementation 
in  complex  geometries.  Meshes  that  are  predetermined  present  a  significant  constraint. 
Flexible  mesh  distributions  are  easily  extended  to  geometries  that  are  not  tensor  products 
of  straight  lines  (to  be  shown  in  a  later  work). 

Spectral  methods  that  are  not  constrained  to  specific  nodal  points  would  clearly  be 
more  flexible  than  conventional  spectral  methods.  Specifically,  a  distribution  of  points  that 
more  closely  approximates  the  disparate  features  in  the  domain  could  be  adopted  from 
the  outset.  Subsequent  adaptation  to  solution  features  in  the  domain  need  not  rely  on 
smooth  mappings.  In  addition,  these  “arbitrary-grid  spectral  techniques”  could  be  used  in 
conjunction  with  multidomain  ideas.  We  focus  on  formalizing  these  ideas  within  the  context 
of  spectral  techniques. 

In  this  paper,  we  present  some  ideas  for  constructing  spectral  methods  with  arbitrary 
grids.  We  demonstrate  these  ideas  for  a  case  of  spectral  solutions  of  hyperbolic  equations; 
however,  these  ideas  can  be  applied  to  any  partial  differential  equation.  To  illustrate  the 
basic  idea,  consider  the  following  hyperbolic  system  of  equations  in  conservation  form: 

with  arbitrary  initial  and  boundary  conditions.  For  spectral  methods,  a  polynomial  (in  the 
spatial  variable  x)  of  degree  N,  Uiv(x,t),  and  a  projection  operator  /^v  are  sought  such  that 

,  fdUjv  dIjvF([/jv)]  . 

Of  the  spectral  techniques,  the  most  popular  method  is  the  Chebyshev  collocation  method, 
in  which  In  fix)  collocates  fix)  at  the  Chebyshev  GL  points  —  cos(|^).  Note  that  we 


have  here  two  projections;  one  involves  the  differentiation  of  F{Un),  and  the  other  involves 
the  way  that  the  equation  is  satisfied.  Thus,  the  first  application  of  the  operator  In  occurs 
when  we  approximate  by  the  derivative  of  the  interpolation  polynomial  that  interpo¬ 

lates  F{U)  at  the  Chebyshev  GL  quadrature  nodes.  The  second  application  of  In  occurs 
when  we  satisfy  the  approximate  equation 

'dUN  _  dlNFjUN)]  ^  Q 
dt  dx 

at  the  Chebyshev  GL  points. 

The  basic  premise  for  unstructured  spectral  methods  is  that  equation  (2)  does  not  have  to 
be  satisfied  in  the  same  manner  in  which  the  operation  is  carried  out.  In  particular, 

the  derivative  operation  can  be  carried  out  by  interpolation  at  any  selected  points;  the 
equation  is  satisfied  by  either  a  Galerkin  formulation  or  by  a  collocation  method  at  a  different 
set  of  points.  Most  importantly,  the  equation  must  be  satisfied  correctly. 

Mathematically  speaking,  we  can  replace  equation  (2)  with 

,  \dUN  dJNF{UN)]_^ 

[“5( - & — J  =  “ 

where  In  ^  Jn- 

In  reference  [4],  a  particular  case  with  this  approach  has  been  discussed.  The  operator 
Jn  was  defined  by  the  Chebyshev  collocation  operator,  and  In  was  the  Legendre  collocation 
operator.  In  the  constant-coefficient  case  {F(U)  —  U),  this  method  reduces  to  the  Legendre 
collocation  method  with  an  efficient  way  of  calculating  the  derivative  by  using  Chebyshev 
collocation  points.  We  now  generalize  this  notion  to  an  arbitrary  set  of  points,  which  enables 
us  to  apply  spectral  methods  in  circumstances  for  which  the  grid  points  are  not  nodes  of 
some  Gauss  quadrature  formula. 

The  method  discussed  in  this  paper  is  different  from  using  a  transformation  to  redis¬ 
tribute  the  grid  points.  The  use  of  a  transformation  to  redistribute  the  grid  points  involves 
approximation  of  the  solution  by  a  polynomial  in  the  transformed  variable;  as  a  result,  the 
approximation  is  not  a  polynomial  in  the  original  variable.  Our  method  utilizes  a  polynomial 
in  the  original  variable.  Moreover,  the  new  method  can  be  applied  to  totally  unstructured 
grids. 

Finally,  it  should  be  noted  that  the  new  method  has  many  similarities  with  spectral 
elements,  although  the  method  of  derivation  is  different.  For  instance.  Patera  [5]  or  Korczak 
et.  al  [1]  used  global  polynomials  (Lagrangian  interpolants),  passed  through  the  Chebyshev 
collocation  points,  to  obtain  spectral  elements.  However,  their  work  was  not  generalized  to 
arbitrary  grids. 
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2  The  Differentiation  Matrix  for  Unstructured  Grids 


Consider  the  set  of  points  (xq  =  1,X\,X2, xyv-i,  =  ~1))  where  the  points  xi,X2-,  •••,  x^-i 
are  arbitrary.  Let  f{x)  be  a  function  defined  everywhere  in  [—1,1].  The  interpolation 
polynomial  fN{^)  fbat  collocates  f{x)  at  the  points  Xj  is  given  by 

/„(!)  =  j«/  =  X;/(i,)i,(x)  (4) 

j=o 

where  the  Lagrange  polynomials  Lj{x)  are  defined  by 


L{x)  =  {x  -  Xi){x  -  X2)...{x  -  XN-2){x  -  XN-i) 
(1  —  x^)L{x) 

=  {l-x]){x-xj)L'{xj) 

"  2L(1) 

T  (r\  {'i--x)L{x) 


The  Lagrange  polynomial  evaluated  at  the  discrete  points  Xk  for  k  ^  j,  is  equal  to  0; 
Lj{xk)  —  ^j,k  • 

We  use  as  the  approximation  to  Note  that  ^  has  two  alternative  repre¬ 

sentations;  the  first  is  obtained  by  differentiating  (4)  as 

^  =  !;/(*, )i'(x)  (9) 


The  second  representation  stems  from  the  fact  that  is  a  polynomial  of  degree  TV  —  1; 

therefore, 


(10) 

j=0 

Equations  (9)  and  (10)  are  used  to  relate  the  grid-point  values  of  the  derivative  /^{xj) 
to  those  of  the  function.  The  most  obvious  way  is  to  equate  the  expressions  in  (9)  and  (10) 
at  the  grid  points  Xk  {0  <  k  <  N)  to  obtain: 

N 

fNi^k)  =  J2^'ji^k)f{Xj)  (11) 

j=0 

To  rewrite  expression  (11)  in  matrix  form,  we  first  denote 

/'  =  IfNi^o),  fN{xN)f ,  /=  [flxo),...,f{xN)f 
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which  yields 


/'  =  Df  (12) 

where  the  differentiation  matrix  D  is  given  by 

D  =  {dj,k)  =  [L'k{x^)]  (13) 

Another  method  for  expressing  the  equivalency  between  (9)  and  (10)  is  to  state  that  the 
difference  between  these  expressions  (which  is  identically  0)  is  orthogonal  to  all  polynomials 
of  degree  <  N: 


Lk{x)dx  —  0 


0  <  ^  <  AT 


The  system  of  equations  that  follows  from  (14)  can  be  rewritten  as 

N  ^  N 

Z  0<k<N 

j=0  j=0 


where 


mk,j  =  Lj{x)Lk{x)d2 

Sk,j  =  J  ^L'j{x)Lk{x)dx 


In  the  matrix  form,  equation  (15)  becomes 


vhere 


Mf  =  Sf 


M  =  (mkj)  0<j,k<N 


S  —  i^k,j)  0  j,  k  N 

Equations  (14)  and  (11)  are  different  manifestations  of  the  same  fact:  (9)  and  (10) 
equivalent.  Therefore,  the  differentiation  matrices  derived  from  (14)  must  be  the  sam 
the  matrix  derived  from  (11)  (with  the  assumption  that  M  is  invertible): 


same  as 


To  prove  this  directly,  we  show  that 


MD  =  S  (22) 

By  writing  the  (i,  k)  term  on  the  left-hand  side  of  (22),  we  obtain 

N 

^  rrii^jdj^k 
j=o 

If  we  substitute  (13)  and  (16)  into  (2),  then  we  get 

(MT>).,=  /'  Li{x) 

•'-1  [j=:0 

We  now  use  the  fact  that  every  polynomial  of  degree  N  is  identical  with  its  N -degree  inter¬ 
polation  polynomial.  Thus,  because  L'i^{x)  is  a  polynomial  of  degree  TV  —  1  and 

Y,L[{x,)Lj{x) 

j=o 

is  its  interpolant  at  the  points  Xj  {0  <  j  <  N)  then 

N 

Y,Lji^)L'kixj)  =  L[{x) 

j=0 

which  yields 

(MB)j  ^  =  £  ii(x)ll(x)  =  Si,t 
(which  is  apparent  from  (17)).  This  establishes  expression  (22). 

□ 

Thus,  we  have  defined  a  new  method,  based  on  the  arbitrary  distribution  of  points,  to 
approximate  the  derivative  of  a  function.  The  attractive  features  of  the  representation  (21) 
of  the  differentiation  matrix  are  summarized  in  lemma  3.1  and  lemma  3.2; 

Lemma  3.1: 

The  matrix  M  defined  in  (16)  is  a  symmetric  positive- definite  matrix. 


5 


Proof: 


The  fact  that  M  is  symmetric  follows  immediately  from  the  definition  (16).  In  fact, 

mk,j  =  Lj{x)Lk{x)dx  =  rrij^k 

We  must  show  that  M  is  positive  definite.  Let  V  be  an  TV  +  1  component  vector: 


V  =  {Vo,...,Vn) 


Then, 


N  N 

v^Mv  = 

i=0  j—0 

Recall  the  definition  of  rrii^j  from  (16).  We  get 


rriijViVj 


MP  =  /  X]  ^  VjLj{x)dx  > 

-'1  -•  rv 


Clearly,  the  equality  sign  holds  only  if  V  is  the  null  vector. 


Equation  (23)  can  be  interpreted  in  a  different  way.  Let  u(a:)  be  the  polynomial  of  degree 
TV  defined  by 


so  that 


Then,  (23)  can  be  rewritten  as 


v{xj)  =  Vj  0<j<N 


^{x)  =  Y.VjLj{x) 


V^MV 


Thus,  every  vector  V  can  be  identified  with  a  polynomial  v(a;)  that  takes  the  values  of  its 
components  at  the  grid  points  Xj.  The  vector  space  norm 


V^MV 


is  equivalent  to  the  function  space  norm 


l\{xr 


Next,  we  will  consider  the  properties  of  the  matrix  S. 
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Lemma  3.2; 


Let  S  be  defined  in  (17),  and  let  V  be  defined  as  before.  Then, 

=  \(vi  -  vf,) 


(25) 


Proof: 

We  start  by  showing  that  S  is  almost  antisymmetric.  ^From  the  definition  (17) 

=  J  ^L'j{x)Lk{x)dx 

and  integration  by  parts,  we  get 

Sk,j  =  ■^i(l)-^*(l)  ~  ~ 

We  now  use  the  definition  of  the  Lagrange  polynomials  (6)-(8)  and  note  that  the  bound¬ 
ary  terms  vanish  for  0  ^  j,k  ^  N  to  yield 

^kj  ~t"  ^j,k  —  dk^Ndj,N 

Thus, 


N  N 

r’-sr  =  EE 

k=0  j=0 
N  N 

=  0  E  E  +  sj,k) 


1 


k^0j=0 
N  N 


=  JEE  VjVkihfidjfi  -  h,Ndj,N) 


^  k=0j  =  0 


which  completes  the  proof  of  (25). 


□ 


As  before,  equation  (25)  has  a  natural  interpretation  in  the  polynomial  space.  Let  v{x) 
be  the  polynomial  of  degree  N  such  that  v(xj)  —  vj.  Then, 

N  N 

v^sv  =  EE  VjVkSkJ 

k=0  j=0 
N  N 


=  f  Lk{x)L'j{x)d. 

k=0  j=0 

=  j  v{x)v'{x)dx 


Note  that 


u(l)  =  uo,  v(-l)  =  vn 

Thus,  (25)  is  an  integration-by-parts  formula. 

The  last  issue  that  we  will  discuss  in  this  section  is  the  relationship  between  differentiation 
matrices,  based  on  different  grid-point  distributions.  Consider  two  grids  Xj  and  yj  {j  — 
0,  ...,A’).  Let  the  Lagrange  polynomial  TJ(x)  be  defined  as  in  (6)-(8),  and  let  Tj(x)  be 
defined  in  a  similar  way,  based  on  the  set  of  points  yj.  This  defines  two  differentiation 
matrices  (see  (11)): 

(26) 
(27) 


Dr  =  (dl^  =  [(if)'(Xj) 


and 


D,  =  =  [(i|)'(w) 

We  now  show  that  the  two  matrices  are  similar. 


Theorem  3.1: 


Define  the  matrix  T  by 

II 

II 

(28) 

Then  define 

1 

(29) 

[ 

D,  =  TD^T-' 

(30) 
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Proof: 


1.  Because  is  a  polynomial  of  degree  N, 


i=0 


If  we  substitute  x  =  then  we  get 


^'j{ym)Ll{Xj)  —  L\{ym)  =  6k,r 


which  proves  1. 


2.  Again,  the  Lagrange  polynomials,  based  on  the  grid  points  t/^,  are  polynomials  of  degr 
A/’;  therefore,  their  derivative  can  be  represented  as 


By  the  same  token. 


(3: 

/=0 

Now,  we  substitute  x  =  y^a  in  (31)  and  (32)  to  get 

{^){ym)  =  X]5Z-^f(2^j)(-^i)'(®0-^f(?/m)  (31 

i=o  /=o 

The  left-hand  side  is  the  (m,  i)  element  of  D„  whereas  the  right-hand  side  is  the  (m, ; 
element  of  T-^D^T-,  thus,  (30)  has  been  proved. 


3  The  Legendre  Galerkin  Method  Based  on  Arbi- 
trary  Grids 

Consider  now  the  linear  form  of  (1): 


Ut{x,t)  =  U,{x,t) 
U[x,{i)  =  f{x) 

=  g{t) 


1  <  X  <  1 


We  introduce  a  new  method  for  the  discretization  of  (34),  based  on  the  differentiation 
matrix  introduced  in  the  last  section.  Note  that  the  differentiation  matrix  uses  the  arbitrary 
grid  Xj.  With  the  new  method,  we  seek  a  vector 

that  satisfies 

du 

M—  =  Su  -  Teo[no  -  5'(0]  (37) 

where 

eo  =  (l,0,0,...0)^ 

The  discussion  on  imposing  the  initial  condition  is  deferred  until  later  in  the  paper  because 
of  subtle  issues  that  involve  convergence.  Here,  we  generally  will  not  use 

“i(0)  =  0  <  i  < 

unless  the  grid  points  Xj  have  special  properties. 

The  structure  of  the  matrices  M  and  S,  indicated  in  (34)  and  (37),  leads  immediately  to 
the  following  stability  result: 


□ 


The  stability  result  (39)  can  be  represented  in  a  different  way  in  view  of  the  equivalency 
between  vectors  and  polynomials  established  in  (24).  Specifically,  let  UN{x,t)  be  an  iVth- 
degree  polynomial  such  that 


UN{Xj,t)  =  Uj{t)  0  <  j  < 


Then,  from  (24)  we  see  that 
1  d  /-i 


=  jW-O-Tf? 

=  ~  w;v(-l,t)^]  -  rUAr(l,t)^ 

Thus,  for  the  polynomial  UN{x,t)  we  have  stability  in  the  usual  L2  norm. 

Now,  we  examine  equation  (37)  from  yet  another  point  of  view.  By,  multiplying  (37)  by 
M-^  we  get 

du  ,  , 

—  =  M  ^Su-rM  ^eo[uQ-g{t)]  (41) 

or  in  view  of  (21),  we  obtain 


—  =  Du-tM  g{t)]. 

The  expression  M~^eo  can  be  evaluated  explicitly. 


Theorem  4.2:, 


Let  M  be  the  mass  matrix  defined  in  (16).  Define  the  residual  vector  r  by 

M-^e-‘o  =  r  =  (r„...,r^)2’ 

Then, 

_  Pn-\-i{^3)  +  -Pjv(^l) 

^  2 

where  Pn{x)  is  the  Legendre  polynomial  of  order  N. 
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Proof: 

We  must  verify  that  if  r  satisfies  (43),  then  the  expression 

Mr  =  Co 

is  also  satisfied.  Substituting  (16)  into  expression  (3)  yields 

N 

(Mr).  = 

j=o 


Substituting  expression  (43)  into  (44)  yields 


(Mr);  =  /'  i..(x)f;4(= 


^N+li^j)  +  P wi^j) 


Because  and  are  polynomials  of  degree  <  iV,  they  coincide  with  their  A^'th- 

degree  interpolation  polynomials;  therefore, 

'^Lj{x)  ^  ' Pn+i{^)  +  P’n{^)' 

3=0  i  2  J  L  2  . 


SO  that 


(Mf).  =  £i,(x) 


P/v+i(^)  +  P 


=  Li{l)  PnI^)  Pjv+i(-l)  +  Pjv(-l) 

2  2 


/>:w 


PAr+i(a:)  +  Pn{x) 


Recall  that 


PMl)  =  l,Piv(-l)  =  (-l)'^ 

and  that  Pjv  and  Pv+i  are  orthogonal  to  all  polynomials  of  degree  <  iV;  the  last  two  terms 
in  the  right-hand  side  of  (46)  vanish,  and  we  are  left  with 

(Mf)^.  =  L,{1)  =  6,,o 


which  proves  theorem  4.2. 
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Theorem  4.2  sheds  a  new  light  on  the  connection  between  the  method  defined  in  (37) 
that  uses  the  arbitrary  set  of  grid  points  Xj  and  the  Legendre  Galerkin  method.  They  are 
the  same  method. 

Theorem  4.3: 

The  method  defined  in  (37)  is  equivalent  to  the  Legendre  Galerkin  method. 


Proof: 


Define 


N 


UN{x,t)  = 

3=0 


where  Uj{t)  are  the  elements  of  u  defined  in  (37).  Then,  UN{x,t)  satisfies  the  error  equation 


5ujv(x,t)  duN{x:t) 


dt 


dx 


—  T 


+  Pn{^) 


[UN{l,t)  -  git)] 


(46) 


The  error  equation  is  satisfied  because  both  sides  of  expression  (46)  are  polynomials  of  degree 
N.  The  two  sides  agree  at  +  1  points  Xj  {j  =  0,  ...,iV)  by  virtue  of  (37),  which  indicates 
that  they  are  equivalent.  Because  the  right-hand  side  is  orthogonal  to  all  polynomials  of 
degree  N  that  vanish  on  the  boundary  x  =  1,  this  error  equation  is  the  same  equation  that 
is  satisfied  by  the  Legendre  Galerkin  method  [6]. 


□ 


As  equation  (46)  demonstrates,  the  precise  method  for  imposing  the  boundary  conditions 
affects  the  overall  behavior  of  the  method.  Section  3  shows  that  two  differentiation  operators 
defined  on  different  grids  are  similar  and,  thus,  have  the  same  eigenvalues.  We  now  show 
that  the  modified  differentiation  matrix  also  has  this  property.  Equation  (42)  produces  a 
modified  differentiation  matrix  (i.e.,  a  differentiation  matrix  that  takes  into  account  the 
boundary  conditions): 

D-tR 

where  the  boundary  matrix  R  is  defined  as 

Ri,j  =  TiSjfi  (47) 
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Suppose  now  that  we  have  two  grids  Xj,  yj  {j  =  0,...,N).  We  have  shown  in  theorem  2.3 
that  Dx  and  Dy  are  similar; 

Dy  =  TDxT-'^ 

where  the  matrices  T  and  are  defined  in  (28)  and  (29).  We  show  now  that  the  same 
similarity  transformation  exists  for  the  modified  differentiation  matrices.  That  is, 


Dy-TRy  =  T{Dx-TRx)T- 


(with  theorem  3.1) 


Ry  =  T  RxT' 


Consider  element  {i,j)  of  the  right-hand  side: 


iV  iV 


l—O  m=0 


EE  myi)riSo,mL] 

l—Q  m=0 


We  recall  that 


^  2 


Rn{xi) 


where  RnIx)  is  a  polynomial  of  degree  N  and  is,  therefore,  equal  to  its  interpolant.  Thus, 

(TRxT  ^  =  i?7v(j/!)T|(a;o)  =  RN{yi)6j^Q 

which  proves  that  the  similarity  transformation  is  valid  even  for  the  modified  derivative 
matrix. 

The  Legendre  Galerkin  method  defined  by  equation  (37)  is  stable;  therefore,  the  initial 
error  is  not  amplified.  However,  the  effects  of  initial  conditions  must  be  carefully  taken  into 
account.  We  know  that  polynomials  based  on  arbitrary  grid  distributions  may  generally  be 
nonconvergent  (the  Runge  phenomenon). 

The  initial  error  can  be  decreased  with  the  number  of  mesh  points  N  by  constructing  the 
Chebyshev  interpolation  as  an  initial  condition.  Thus,  let 

0<)  <N  (49) 

L^(x)  =  {l-x^)TUx)  (50) 

The  Chebyshev  approximation  for  the  initial  condition  is,  then. 


CnUx)  =  £/(6)Tj(a;) 


(52) 


so  that  the  recommended  initial  approximation  will  be 

/(^i)  ~ 

j=0 

This  approximation  will  provide  a  convergent  approximation  for  the  initial  condition.  Of 
course,  the  Chebyshev  approximation  is  not  the  only  possibility;  any  other  spectral  or  pseu- 
dospectral  approximation  would  do  as  well. 

We  now  briefly  discuss  the  issue  of  implementation.  Two  methods  are  available  for 
implementing  the  arbitrary-grid  spectral  methods.  The  first  method  is  to  form  the  matrices 
M  and  S  by  carrying  out  explicitly  the  integrations  in  (16)  and  (17).  (This  technique  is 
utilized  in  the  two  examples  presented  later  in  the  text.)  This  procedure  is  done  once  and 
for  all  for  every  given  set  of  grid  points.  Then,  the  equations  are  solved  as  described  in 
(37).  A  mor-e  convenient  method  that  does  not  involve  evaluating  integrals  is  to  use  the 
differentiation  matrix  D  defined  in  (13)  and  solve  the  system  (42)  with  the  identity 

A/T-l;?  _  -^iV+lC^j)  + 

M  eo-  ^ 

proven  in  theorem  4.2.  For  a  large  N,  the  method  that  will  be  the  most  successful  is  the  one 
with  the  least  sensitivity  to  round-off  errors.  This  point  has  not  been  fully  investigated  at 
this  time. 

Finally,  an  observation  in  regard  to  the  maximum  allowable  time  step  for  the  arbitrary- 
grid  spectral  schemes.  All  spatial  operators  have  the  same  eigenvalues,  regardless  of  the 
spatial  distribution  of  points  (48).  Therefore,  the  maximum  allowable  time  step  is  the 
same  for  all  schemes.  Stability  is  a  matrix  property,  and  depends  on  all  the  points  in 
the  distribution.  This  observation  is  somewhat  counter  to  the  conventional  finite-difference 
notion,  in  which  the  maximum  time  step  is  governed  by  the  smallest  grid  spacing. 

4  The  Legendre  Collocation  for  Unstructured  Grids 

The  Legendre  collocation  for  unstructured  grids  involves  the  approximation  of  the  integrals 
in  (16)  and  (17)  by  the  Gauss-Lobatto-Legendre  (GLL)  quadrature  formula.  Let  {tjq  = 
l,?7i,...,r/iv-i,??Ar  =  —1)  be  the  nodes  of  the  GLL  quadrature  formula  and  loi,  0  <  I  <  N  he 
the  weights.  We  define  a  new  mass  matrix  Me  by 

JV 

Mc(f,j)  =  (53) 

/=0 

where  the  Lj{x)  are  the  Lagrange  polynomials  at  the  points  (xq  =  T  a^i,  a;2, ...,  x^r  = 
— 1).  Note  that  this  is  an  arbitrary  set  of  grid  points. 
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The  matrix  may  be  different  from  M  because  the  GLL  formula  is  exact  to  order  2N—1 
and  Lj{x)Lk{x)  is  a  polynomial  of  order  2N.  The  matrix  M,;  is,  however,  a  symmetric  and 
positive-definite  matrix. 

By  introducing  quadrature  to  equation  (17),  we  define  a  new  stiffness  matrix  as 

N 

(54) 

/=0 

Note  that  because  of  the  exactness  of  the  GLL  formula,  the  sum  on  the  right-hand  side  of 
(54)  is  the  same  as  the  integral  in  the  right-hand  side  of  (17);  therefore, 

Sc  =  S 

For  this  reason,  the  property  (25)  is  true  for  the  stiffness  matrix  Sc  also. 

The  uniqueness  of  the  differentiation  matrix  D  also  yields 

m;'Sc  =  M-^s 

which  does  not  contradict  the  fact  that 

Mc/M 

because  the  matrices  Sc  and  S  are  singular. 

In  the  Legendre  collocation  method  of  (34)  with  arbitrary  grids,  we  seek  a  vector 

u  =  [uo(t),...,U7v(t)]^ 

that  satisfies 


Me—  =  ScU  -  Teo[uQ  -  g{t)]  (55) 

where 

eo  =  (1,0,0,...0)^ 

Alternatively, 

du 

—  =  Du-TM;^eo[uo-g{t)]  (56) 

The  stability  of  (55)  follows  immediately  from  the  fact  that  Me  is  symmetric  positive 
definite  and  Sc  satisfies  (25).  Our  aim  is  to  show  that  (55)  is  equivalent  to  the  usual  Legendre 
collocation  method. 


16 


Theorem  4.1: 


Let  Me  be  the  mass  matrix  defined  in  (53).  We  define  the  residual  vector  r  by 


Then, 


M/eo  =  f  =  (re,...,rAr)^ 


where  Pn{x)  is  the  Legendre  polynomial  of  order  N. 


Proof: 


We  start  by  noting  that  the  nodes  t}i  of  the  GLL  formula  are  the  zeroes  of  the  polynomial 

P;(a:)(l  - 

Also,  because  Pjv(l  +  x)  is  a  polynomial  of  degree  iV, 

j=0 


Therefore, 


i=o 
N  N 

=  E  E  LjiVl)Li{Vi)PN{xj)il  +  Xj)u;i 

j=0  /=0 
N 

=  Yl^iiVl)PNiVl){'^  +  Vl)<^l 


=  2N^6o,; 


which  proves  the  theorem. 


Equation  (56)  also  can  be  viewed  as  shown  in  the  following  example.  We  seek  a  polyno¬ 
mial  (in  x)  UN{x,t)  of  the  form 


Ujv(a;,<)  =  'Y^UN{xj,t)Lj{x) 

i=o 
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(58) 


such  that 

— =  '^UN{xj,t)L'j{xk)  -  TRN{xk)[uN{l,t)  -  git)]  (59) 

j=o 

where 

Rjv(x)  =  (1  +  x)F^(x) 

This  approach  is  equivalent  to  the  Legendre  Collocation  Method  [6]. 

The  extension  of  the  arbitrary-grid  Legendre  collocation  method  from  the  linear  case  (34) 
to  the  solution  of  the  nonlinear  case  (1)  is  immediate.  The  issue  of  implementation  could  be 
significant.  To  avoid  computing  the  points  qi,  the  best  choice  is  to  use  the  formulation  (56) 
rather  than  (55).  In  this  case,  Me  and  Sc  do  not  need  to  be  computed. 

At  this  stage,  note  that  for  the  case 

RNix)  =  Pn{x) 

we  have  the  Legendre  Tau  method,  with  the  additional  property  of  an  improved  time  step. 
However,  we  do  not  have  the  representation  of  the  Legendre  Tau  method  in  the  form  of  (55). 


5  Unstructured  Grids  for  Unbounded  Domains:  La- 
guerre  Methods 


Consider  the  equation 


dt 

UiO,t) 

Uix,0) 


dx 

9ii) 

hix) 


0  <  X  <  oo 


Note  that  the  domain  is  semibounded.  Note  also  that  if  git)  =  0,  then 

—  /  e-^U\x,t)dx  =  -  /  e-^U^ix,t)dx 
at  Jo  Jo 

Assume  that  we  have  an  arbitrary  set  of  grid  points 


(60) 

(61) 


(62) 


(Xo  =  0,Xi,...,XAr) 


In  the  Galerkin  procedure,  we  approximate  the  derivative  of  a  function  fix)  whose  values 
at  Xj  are  given  by  the  derivative  of  its  interpolant  fNix).  After  we  define 

T(x)  =  (x  —  xo)...(x  —  xn) 
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(63) 


Equation  (67)  defines  the  differentiation  matrix  D.  In  fact,  if  we  define 

rufcj  =  (^Lj,Lk)  (68) 

and 

(69) 

where  the  scalar  product  {u,v)  is  defined  as 

1*00 

(u,  u)=  /  e~^u{x)v{x)dx 
Jo 

then  we  get 

D  =  M-^S 

As  before,  the  differentiation  matrix  is  unique.  The  manner  in  which  the  matrices  M  and  S 
are  constructed  leads  immediately  to  the  following  lemma. 

Lemma  6.1: 

The  matrix  M  is  symmetric  positive  definite.  The  matrix  S  satisfies 

S  +  =  M  —  diagonal(l,0, 0,...,0)  (70) 
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Proof: 


i.From  the  definition  of  the  matrix  S,  we  have 

^k,j  —  {Lj^Lk) 

=  /  e  ^ L  ■{x)Lk{x)dx 

Jo 

roo 

—  —Lj{0)Lk{0)  —  e  ^ Lj{x)Li^{x)dx 

J  0 

roo 

+  /  e~^Lj{x)Lk{x)dx 

Jo 

By  using  the  definition  of  the  matrix  M  and  the  properties  of  the  Lagrange  polynomials,  we 
get 

^k,j  ^j,k  “t"  ^k,j 

which  proves  (70). 

□ 

To  discretize  (60),  we  introduce  the  unknown  vector 

u  =  [uo{t),...,UN{t)f 

that  satisfies 

_  _  dxi  „  ^  \  T 

M—  =  -Su  -  Teo[t<o  -  5^(0]  (72) 

The  stability  is  immediate,  as  shown  in  the  following  lemma. 

Lemma  6.2: 


Let  u  satisfy  (72),  with  g{t)  =  0.  Then,  we  have  the  energy  estimate 


-■if'Mu  =  —xFMu  —  (2r  —  1)^^ 


Proof: 

Equation  (73)  follows  immediately  from  multiplying  (72)  by  and  using  (70). 


□ 
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Lemma  6.1  implies  that  the  method  is  stable,  provided  that  t  >  |.  Note  that  the  energy 
estimate  (73)  for  the  approximation  is  nearly  the  same  as  for  the  differential  equation  (64). 

We  still  must  show  that  the  method  described  in  (72)  is  equivalent  to  the  Laguerre 
Galerkin  method.  We  begin  by  rewriting  (72)  as 

^  g{t)]  (74) 


The  key  issue  is  to  identify  the  vector 


which  is  done  in  the  following  theorem. 


Theorem  6.1: 


M-'eo 


Let  M  be  the  mass  matrix  defined  in  (68).  Define  the  residual  vector  f  by 


Then 


M  ^eo  =  r  =  (ro,...,rAr)^ 


_  ±A0)  I 

^  ~  dx 


where  is  the  Laguerre  polynomial  of  order  N. 


Proof: 


We  must  verify  that  r  satisfies  (75)  and  that 


Mr  =  eo 


We  begin  by  expanding  (Mr)  as 


(Mr),  =  E 


yoo  _  N 

=  e  ^Li{x)  E  Lj{x)rjdx 

•'O  7=0 


If  we  substitute  (75)  into  (76  ),  then  we  get 


(Mr),  =  r 
Jo 


j=o  L 


Because  is  a  polynomial  of  order  N,  it  coincides  with  its  interpolant;  therefore, 


Thus, 


(Mf);  =  /' 

Jo 


e  "" Li{x)—C%\-^{x)dx 


If  we  integrate  the  right-hand  side  by  parts,  we  get 


/’oo  r  oo 

(Mr).  =  -A(0)£);;^i(0)  /  e-^Li{x)6^\^{x)dx  -  /  L\{x)C%\,{x)dx 

JO  Jo 

The  last  two  terms  on  the  right  vanish  because  of  the  orthogonality  of  and  the  first 
term  vanishes  if  z  ^  0;  thus, 

(Mr-)^.  =  -<5,,o 

and  the  theorem  is  proven. 


Another  method  for  getting  the  Lageurre  method  on  the  grid  Xj  is  to  seek  a  polynomial 
UN{x,t)  such  that 


duNiXk,t)  _  ^l’.(^Xk)uN{Xj,t)  -T£%{xk)[uN{xo,t)  -  g{t)]  (78) 

i=o 

where  {LG)Nix)  is  the  Mh-degree  Laguerre  polynomial.  This  approach  is  the  Laguerre 
collocation  method. 


6  Numerical  Results 


We  now  test  the  previous  theoretical  results  with  two  numerical  examples.  The  linear  equa¬ 
tions  (.34)-(36)  are  solved  with  f{x)  =  sin(7ra:),  g{t)  =  sin[7r(l  -1- 1)],  and  the  exact  solution 
U{x,t)  =  sin[7r(a;-}-t)].  A  variety  of  grids,  from  Chebyshev  to  “randomly  generated”  grids, 
are  used  to  test  the  accuracy  and  stability  of  the  method.  For  all  calculations,  128-bit 
arithmetic  is  used  to  ensure  adequate  precision. 
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Figure  1  shows  the  refinement  study  on  five  different  grids: 


1.  Uniform  grid  Xj  =  {j  =  0,...,iV) 

2.  Chebyshev  grid  Xj  =  cos{^) 

3.  A  linear  combination  of  the  uniform  grid  and  Chebyshev  grid  (i.e,  Xj  =  + 

0.5cos(f)) 

4.  Chebyshev^  (i-^.,  Xj  =  cos^(^)) 

5.  (Chebyshev)^  for  — 1  <  a;  <  0  and  (Chebyshev)  2  for  0  <  x  <1,  where 

(Chebyshev)2  fo  defined  by  the  grid  points  xj  =  cos2(^). 

The  logio  of  the  L2  error,  plotted  against  the  number  of  points  in  the  approximating 
polynomials  is  shown  in  Figure  1.  The  problems  are  run  to  the  physical  time  T  —  2.  The 
convergence  is  exponential  for  all  cases  until  machine  round  off  is  encountered.  These  results 
are  consistent  with  the  previous  numerical  results.  (Note  that  the  Chebyshev  grid  is  the 
least  sensitive  to  round  off.) 

The  Legendre  Galerkin  method  defined  by  equation  (37)  is  stable;  therefore,  the  initial 
error  is  not  amplified.  However,  the  effects  of  initial  conditions  must  be  carefully  taken 
into  account.  We  know  that  polynomials  based  on  arbitrary  grid  distributions  generally 
may  be  nonconvergent.  This  property,  called  the  Runge  phenomena,  is  easily  demonstrated 
by  approximating  the  function  f{x)  =  ^  <  1)  on  a  uniform  grid.  The 

global  approximating  polynomials  oscillate  wildly  at  each  end  of  the  domain,  which  yields 
a  poor  approximation  in  those  regions.  The  Runge  phenomena  is  alleviated  by  using  a  grid 
distribution  (like  the  Chebyshev  grid  distribution),  which  clusters  points  near  the  boundaries 
—  1  and  1. 

Figure  2  illustrates  that  a  Runge-like  phenomena  exists  within  the  arbitrary-grid  spectral 
methods  if  special  precautions  are  not  taken  in  the  initialization  step.  In  this  problem,  the 
linear  equations  (34)-(36)  are  solved  with  f{x)  =  g{t)  =  :j7j]2 ;  exact 

solution  t/(a:,  t)  =  •  The  simulation  is  run  to  time  T  =  0.001  (a  physical  time  that 

occurs  well  before  the  influence  of  the  initialization  is  lost.)  (Running  to  a  physical  time 
T  >  2  yields  exponential  convergence  on  all  grids.)  Convergence  is  achieved  only  for  the 
Chebyshev  grid  distribution. 

The  source  of  the  error  in  this  problem  is  the  failure  of  the  arbitrary  grid  that  approx¬ 
imates  the  polynomial  to  converge  to  the  initial  condition.  For  small  times  (less  than  1 
convective  sweep),  erroneous  information  is  left  in  the  domain,  and  the  resulting  method  is 
nonconvergent.  By  changing  the  problem  slightly,  however,  convergence  can  be  recovered  on 
all  grids. 
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To  initialize  the  problem,  we  must  construct  an  approximation  to  the  initial  condition 
/(x),  based  on  the  grid  points  Xj  (0  <  j  <  N).  We  want  to  keep  the  flexibility  and  rigid 
structure  of  the  original  grid  distribution;  however,  the  interpolation  polynomial,  based  on 
the  grid  points  Xj,  generally  is  not  convergent.  Therefore,  we  use  the  method  outlined  in 
(49)  and  (52).  With  this  initialization,  spectral  convergence  is  recovered. 

7  Conclusions 

A  new  technique  for  implementing  spectral  methods  for  hyperbolic  equations  has  been  devel¬ 
oped  that  does  not  require  grid  points  that  are  nodes  of  some  Gauss  quadrature  formula.  For 
this  reason,  this  method  is  referred  to  as  an  arbitrary-grid  spectral  method.  Both  Galerkin 
and  collocation  formulations  are  presented  for  the  conventional  Legendre  method,  and  a 
Galerkin  formulation  is  presented  for  the  conventional  Laguerre  method. 

The  basis  for  the  stability  of  the  unstructured  spectral  schemes  relies  on  a  weighted 
energy  norm  in  all  cases.  Stability  is  proven  for  the  constant  coefficient  hyperbolic  system. 
All  unstructured  spectral  methods  utilize  a  “weak”  imposition  of  the  boundary  condition, 
similar  to  the  technique  used  in  the  penalty  formulations  of  the  finite  element  method.  With 
this  imposition,  the  complete  differentiation  matrix,  including  boundary  conditions,  is  similar 
to  (i.e.,  it  has  the  same  eigenvalues)  the  conventional  differentiation  operator;  therefore,  this 
matrix  behaves  similarly. 

The  new  formulations  are  demonstrated  on  two  scalar  hyperbolic  problems.  The  arbitrary- 
grid  Legendre  Galerkin  method  is  used  in  both  cases.  Exponential  accuracy  is  shown  in  both 
cases  on  arbitrary  grids.  Care  must  be  exercised  in  the  initialization  procedure  to  ensure 
convergence  of  the  new  schemes. 
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FIGURE  2.  Divergence  of  the  arbitrary  grid  Legendre  Galerkin  method  for  improperly  imposed 
initial  conditions. 
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